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Abstract 

We present an Initial Value Representation for the semiclassical coherent state 
propagator based on complex trajectories. We map the complex phase space into 
a real phase space with twice as many dimensions and use a simple procedure to 
automatically eliminate non-contributing trajectories. The resulting semiclassical 
formulas do not show divergences due to caustics and provide accurate results. 
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1 Introduction 



The investigation of wavepackets dynamics by semiclassical methods has practical 
importance for calculations of several processes involving atoms and molecules. It is 
also a fundamental topic in the study of the classical-quantum connection, especially 
for chaotic systems and for open systems, coupled to environments. 

The history of semiclassical methods goes back to the origins of quantum mechan- 
ics itself. One fundamental result is the so called Van Vleck approximation to the 
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coordinate propagator, derived in 1928 [I], that can be written as 
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(1) 



In this expression S(xi,Xf,T) is the classical action of a trajectory connecting coor- 
dinates Xi to Xf in the time T, m qp is an element of the tangent matrix, that controls 
the motion in the vicinity of this trajectory, and k is the number of focal points 
(where m qp goes to zero) along the trajectory. If more than one trajectory satisfying 
these boundary conditions exists, one has to sum their contributions. From this basic 
propagator one can compute the time evolution of arbitrary wavef unctions. 

A more direct approach to calculate the time evolution of wavepackets in given by the 
propagator in the coherent state representation. The coherent states of the harmonic 
oscillator are minimum uncertainty wavepackets and define a representation involving 
both the coordinates and the momenta that can be readily visualized in the phase 
space. The coherent state propagator (zf\e~% HT \zo) represents the amplitude prob- 
ability that the initial wavepacket \z ) centered on go,Po is found at the state \zf), 
centered on qf,Pf, after a time T. However, the direct evaluation of the semiclassical 
limit of this propagator results in an expression bearing the same difficulties of the 
Van Vleck formula [2,3, 4, 5j, namely: (a) the classical trajectories needed are defined 
by mixed initial-final boundary conditions, rendering the calculation hard, specially 
in multidimensional or chaotic systems and; (b) the formula diverges at phase space 
focal points. Moreover, the trajectories are complex and some of them, even satisfying 
the appropriate boundary conditions, lead to unphysical contributions and must be 

Several methods have been developed to overcome these difficulties, most of them 
based on the idea of initial value representations (IVR) [I1^16M17|[T8 , 19,20 21 22,2 3124"] . 
Among these, the Herman-Kluk propagator (18] and the method of linearized cellu- 
lar dynamics developed by Heller and Tomsovic [25|26] stands out as very accurate. 
More recent derivations and corrections to the basic Herman-Kluk formula have also 
provided new insight into this class of approximation [2Tll28f29|l3"0~] . 




In spite of the many difficulties involved in the calculation of the coherent state prop- 
agator with complex trajectories, this approximation turns out to be very accurate 
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[8ll9lll0lll2|[T3] . Recent work on Bohmian mechanics have also employed complex tra- 
jectories, providing a new formulation leading to accurate results [3T1I32] . In this paper 
we propose the construction of an initial value representation for this approximation 
that removes most of its problems: the mixed conditions defining the trajectories are 
replaced by initial conditions; the complex trajectories are mapped into real trajecto- 
ries of an associated Hamiltonian; the divergences due to focal points are eliminated 
and; a simple and automatic filtering is used to eliminate the non-contributing tra- 
jectories. 

This paper is organized as follow. The next section reviews the semiclassical coherent 
state propagator and its semiclassical approximation in terms of complex trajecto- 
ries. In section |3] we develop the initial value representation for complex trajectories 
(CIVR). It presents what we call the sudden CIVR, where only trajectories satisfying 
the original mixed conditions are considered, and the smooth CIVR where the neigh- 
borhood of the relevant trajectories are considered as well. We also discuss how the 
complex trajectories calculations are performed in terms of real trajectories. This sec- 
tion ends with a discussion of the criteria used for filtering out the non-contributing 
trajectories. In section H] the smooth CIVR is applied to an anharmonic quartic os- 
cillator and some final comments are made. In two appendices useful relations of the 
tangent matrix are derived for both real and complex trajectories. 

2 The semiclassical coherent state propagator 

The coherent state \z) of a harmonic oscillator of mass m and frequency u is defined 
by 

\z) = e -^ l V at |0> (2) 
with |0) the harmonic oscillator ground state and 

In these equations q, p, and a) are operators; q and p are real numbers and z is 
complex. The parameters b = (h/muj) 1 and c = [hmuj)" 1 define the length and 
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(3) 



momentum scales, respectively, and their product is h. 



For a time-independent Hamiltonian operator H, the propagator in the coherent 
states representation is the matrix element of the evolution operator between states 
\z ) and \zf) [55] : 



K(z},z ,T) = (z f \e-n H1 \z ) 



(4) 



The semiclassical evaluation of K(zf, Zq, T) was presented in detail in [HE] - The result 
is given by 



K sc (z* f ,z ,T) = ^ 



d 2 S 



\ Hdz dz*j 



exp{l(5 + /)-i(|^| 2 + |z | 



(5) 



where 



S = S{z*j, zq, t) = J dt' 
o 



ih 



(iiv — vu) — H(u, v, t') 



ih 



- -{u{T)z} + zov(0)) (6) 



is the action and the classical Hamiltonian function is calculated from the Hamilto- 
nian operator as H(u,v) = (v\H\u). The term 

T 

(7) 



, 1 > d 2 H 
I = - -^-^dt 



2 J dudv 
o 



is a correction to the action. The sum over v represents the sum over all contribut- 
ing (complex) classical trajectories satisfying Hamilton's equations with boundary 
conditions 



j_(m +i p(py 

V2\ b c 



1 (q{T) .pgr 
v^l b c j 



In all these expressions the variables u and v are defined by 



u 



V2\b 1 c 



V2\b % c 



(9) 



They are manifestly independent (u ^ v* since q and p are complex), and replace z 
and z* to avoid confusion. In these variables the boundary conditions become 



u(0) = z , v(T) = z) . 



(10) 
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3 A complex initial value representation 



3.1 Basic idea 



The first of the boundary conditions (TTUj) specifying the complex trajectory can be 
written explicitly as 

m +Q m.*+i£z, (id 

b n b n 

where go and po define the initial coherent state \zq). This condition is not sufficient 
to determine the trajectory, since g(0) and p(0) are complex. The missing condition 
is given by the second equation in (JTOj) and refers to the final propagation time T. 



In order to avoid dealing with mixed initial-final conditions, let us first suppose we 
have had a second equation of the form 

^-^ = T-«TT- (12) 
b n b n 

By solving for g(0) and p(0) one finds 
<?(0) = \ [(go + qi) + i^(po - Pi 



(13) 



p(0) = i(fpo + pi) + i&(q x - q 



For q and po fixed, each q x and p\ defines a trajectory with end points q(T) and 
p(T). 

Let q± and pi be the values of q± and p± such that the second of equations ( flOl) is 
satisfied, i.e., for which the initial conditions ffT3l) leads to 



SSl. ib ?Sl = SL- i 0. ) (14) 

b n b n 



where qf and pf define the final coherent state \z/). Then, we can rewrite the semi- 
classical propagator of Eq. (JHJ) as 



K sc (z* f , Zo, T)= J dqxdpxdaiq! - qi)5 a (Pi - Pi) 



X 



i d 2 S i(5+/) _i(|, / | 2+ko |2) 



^ hdz dz*f 



e *ls+'J-3VW *+l«ol ^ (15) 



where the trajectories are now calculated according to the initial conditions (fT3j) and 
their contributions filtered out by the delta functions. Since these are sharped peaked 
functions, equation (TT5|) is identical to ([SD, because only the trajectories satisfying 
the proper boundary conditions (TTUj) are taken into account. If the delta functions are 
replaced by Gaussian functions of width a, Fillinov type expansions become possible 
and a smoothed and better behaved expression arises. 

The equivalent expression of the semiclassical propagation for an arbitrary initial 
state described by the wave-function i/j(zq,0) = {zq\i/)) is 



if)(Zf, T) = ^—^ I K sc (z* f , z , T)4>(zq, 0)6 a (q 1 - q\)5 a (p 1 - pi) dq dp dqidpi. (16) 



We note that the second derivative of the action with respect to its arguments, as 
appearing in the pre-factor of the semiclassical propagator, can be written in terms 
of the tangent matrix, that controls the classical motion in the vicinity of a given 
trajectory. In appendix |A] we derive several useful relations between the tangent 
matrix in u, v and q, p variables for complex and real trajectories. In particular, we 
show that 

' 923 1 (17) 



h dz dz} M, 



Before we end this subsection we define the scaled coordinates and momenta q = q/b 
and p = pb/h. Defining the scaled Hamiltonian 

H(q } p) = ^H(bq,hp/b) (18) 



it is easy to check that the semiclassical expressions in terms of q, p and H become 
identical to the original expressions with b and h replaced by 1. Therefore, from now 
one we shall use these scaled variables, which amounts to set b — H = 1, but will omit 
the bar to make the notation simpler. The original variables will be recovered later 
in the examples. 



3.2 The calculation of complex trajectories 



For analytic Hamiltonian functions H(q,p) it is possible to rewrite the equations of 
motion for the complex variables q and p in terms of real trajectories of an auxiliary 
Hamiltonian system with twice as many degrees of freedom, or as we call it, the 
double phase space. The definitions [Mf35] 



q = Qi + iP 2 , 



P = Px+ iQ 2 



(19) 



and 



H(q,p) = H 1 (Q 1 , Q 2 , P u P 2 ) + iH 2 {g x , Q 2 , P t , P 2 ), (20) 

where H\ and H 2 are real functions, allows to show easily that Hamilton's equations 
for q and p are equivalent to 

*=- w,' (21) 

Note that H 2 is also a constant of the motion. The separation of variables in (fT9l) 
may look unusual because it mixes q's and p's, but this is the proper combination to 
get the correct signs in Hamilton's equations. These separation of variables also look 
natural when the form of equation (|T3|) is considered. 

For the case \ip(0)) = \zq), the real trajectory starting from the center of the wavepacket 
plays an important role, and we shall use it as a reference. Therefore the integration 
over qi and p\ in the CIVR will be centered on q and po and only a limited region 
around this point is expected to significantly contribute to the propagation. In this 
way we write 

qi = q + Aq, Pi = Po + Ap (22) 

and the initial conditions (IT51) reduce to 

Qi(0) =q + Aq/2 Q 2 (0) = Aq/2 

(23) 

P 1 (0) =p + Ap/2 P 2 (0) = -Ap/2. 
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In accordance with Eq. (|T9|) . g(0) = qo + w, p(0) = po + iw with w = (Ag — iAp)/2, 
which is exactly the variable used in a search procedure developed in ref. [9J. 

All the tangent matrix elements appearing in equations (155j) and (j2HD can be readily 
computed from the tangent matrix of the real trajectory in the double phase space. 
This procedure eliminates the need to work with complex trajectories and also the 
so called root search problem, involved in finding trajectories with mixed initial- final 
conditions. 

3.3 The connection between initial and final displacements 

The connection between the initial and final displacements can be established as 
follow. Initially by comparing equations (1T91) with (1T31) we see that 

Qi(0) =i(q Q + qi) 

Q 2 (0) = !(<&- g„) 

(24) 

Pi(0) =f(p +Pi) 
P 2 (0) =±(p -;Pi) 

which also leads to gi = Qi(0) + ^2(0) and p\ = -Pi(O) — P^ti). It turns out to be 
convenient to extend this definition to 

qx(t) = Q 1 (t)+Q 2 (t) 

(25) 

p 1 (t) = P 1 (t)-P 2 (t). 

Because of the filtering functions in ffT5]) and f)16p (smoothed or sharp) the relevant 
contributions to the integrals over qi and p\ come from the vicinities of q~\ and pi, 
that should also be close to go and pq. For this particular trajectory v (T) = z*f. 

[Q X {T) + iP 2 (T)} - tiP^T) + iQ 2 {T)} =q f - iq f (26) 
or, according to (1231) . qi(T) = qj and Pi(T) = pf. 



For neighboring trajectories we may expand the final values of qi{T) and pi(T) around 
qj and Pf as: 



<H(T) *Qf + ^(ft - Qi) + ^(Pi - Pi) 

pi eo « p/ + - + ^?(pi - Pi) 



or 



/ \ 

qi(T)-q f 



P\{T)~Pf 



( 



d gi (T) d qi (T) 



\ 



dqi 
dp x {T) 



dpi 
dp x {T) 



\ dqi dpi J 



( \ 

qi - qi 



. Pl "Ply 



A 



/ \ 

qi - qi 

. Pl "Ply 



(27) 



It follows that 



%l-gi)5(pi - Pi) = |detA| 8{q 1 {T)-q } )5{p l {T)-p } ) 



(2f 



In appendix IB1 equation (IB. 6|) . we show that det A = |M, 



5.^ Sudden complex initial value representation 



In the case of sharp delta functions we can use equation 
of our formulas, that we term sudden CIVR. Since 



to write down the first 



V2 



{q l {T)-ip l {T)) = v{T) 



(29) 



and defining 



5 2 (v(T) - z}) = 2v5{ qi (T) - q f )5( Pl (T)-p f ) 



(30) 



we obtain 



r(>(z* f ,T) = [ |M w | 3 / 2 e < ( s+J 5-3(^l 2+ l a °l a )^V'(2S J 0)(J 2 (v(r)-«;) 



crz d 2 vi 



IT IT 



(31) 



where £ is the phase of M vv . Each pair of phase space points qo, po and qi, pi define 
a complex trajectory with initial conditions 

9( n ) = i(?o + gi) + i|(po - Pi) 

(32) 

p(0) = |(p +Pi) +*s(?o 

The contribution of these trajectories to the final result are filtered by the delta 
function. The integration measures are defined as usual as d 2 Zo/ir = dqodpo/27r and 
d 2 Vi/ir = dqidpi/2n. 

Notice that the arguments of S and / in fl3T]) . which were originally (z*f,z ,T), can 
be replaced by (v(T),zq,T), so that both S and / are computed for the trajectories 
defined by (|32|) . Also important is the fact that M vv in the pre-factor has moved from 
the denominator to the nominator so that divergences at caustics are replaced by 
non-contributing trajectories. This is a well known property of IVR's constructed in 
this way. 



3.5 Smooth complex initial value representation 



If the delta functions in the CIVR are replaced by Gaussian functions a more well 
behaved approximation is obtained. Following Filinov [36] and Makri [37] we replace 
the filtering integrals of trajectories according to 



7r J 7ra 2 

Kqi(T)-q f ) 2 +(pi(T)~p f f]d 2 Vi 



g 2a 2 \M vv \' 2 1 



na 2 



M vv 2 e J, 33 

where we have used equation (J27J in the second line and defined the re-scaled width 
a = a\M vv \. (34) 



The use of smooth filters seems appropriate to coherent state propagation. It im- 
plies that not only the trajectories satisfying the exact boundary conditions (fT0|) are 



in 



considered, but also their neighborhood as defined by the parameter a. In this case 
the action S(zj, z , T) in equation (|T5|) cannot be simply replaced by S(v(T), z , T), 
but has to be expanded around each initial value trajectory up to second order. The 
result is 



dS 



1 d 2 S 



S(z* f , z , T) « S(v(T), z , T) + q^W- v(T)) + ^g^yM - v(T)Y 
« S(v(T), z , T) - iu(T)(z} - v(T)) - i^(z} - v(T)) 2 , 



(35) 



where once again we have resorted to expressions derived in appendix iBl 

The smooth CIVR can then be obtained by using equations fl33l) and fl35l) in ffTBT) : 



*P(z},T)= [ \M vv f 2 



exp < 



v(T) - z) 



of 



^(z* ,0) 



d z d v\ 
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(36) 



where 



= i(S + I) + u(T) (z) - v(T)) +^f; (z*f- v(T)) 



M, a , / . ,„,\2 \z f \ 2 \z \ 2 X 



i-. (37) 
2 2 2 y J 



If the initial state to be propagate is itself a coherent state, equation (fToT) . the smooth 
CIVR simplifies to 



K(z* f ,Zo,T) = J |M OT | 3 / 2 exp j< 



v(T) - z) d 2 Vl 



or 



7rcr 



with 



(38) 



= *(S + /) + u(T)z*j + J^L (z) - v(T)) 



2M„ 



,2 \zfl_\zol_X 
2 2 V 



(39) 



In this paper we shall discuss an example of this simple case only 



3.6 Filtering out non- contributing trajectories 

It is well known that not all trajectories satisfying the boundary conditions (fTUl) 
should be included in the semiclassical propagator. The trajectories for which the 

1 1 



real part of the exponent in (1391 is positive must be discarded as they give rise 
to divergent contributions in the semiclassical limit. These trajectories are probably 
associated with forbidden deformations of the integration contours that are necessary 
to derive the semiclassical approximation (j3J). 

For the harmonic oscillator it can be checked explicitly that not only equations ( l38i ) 
and (1391) give exact results but also that Re(<p) < for all complex trajectories. In 
our calculations trajectories satisfying 



where c is a constant, are neglected. We discuss the importance of the cutoff value of 
c in the next section. 

4 Example 

As a simple application of the smooth CIVR we consider the system 



It has been studied also in [TJ] by directly computing the relevant complex trajec- 
tories. The parameters are set to Q = 1, A = 0.4 and h = 1 . For these values the 
ground state energy is E Q ps 0.559 and the first two excited states have E\ ~ 1.770 
and E 2 ~ 3.319. For the initial wavepacket we choose go — 0, Po = —2.0, and b = 1.0. 
This gives E = H(q,p) = 2.0 for the energy of the central trajectory, r m 4.7 for 
its period, and X turn « ±1.6 for its turning points. Figure 1(c) shows a plot of the 
potential function and indicates also the central trajectory energy. 

We momentarily restore the original un-scaled variables to illustrate both the compu- 
tation of the classical Hamiltonian and the scaling process. The classical Hamiltonian 
function is 



Re{4>) > ch, 



(40) 




(41) 



rr 1 2 1 frV2 3A& 2 \ 2 A 4 

H =2 P + 2[ q+ A q + 



( 



n 2 b 2 



4 



+ 



3A6 4 



16 



) 



(42) 
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where b is the width of the wavepacket. In terms of scaled variables (see equation 
({TBI) ) the Hamiltonian becomes 



where to = h/b 2 , v = VL/uj, A = Xh/u 3 and v 2 = v 2 + 3A/2. For the present values we 
have lo = v — 1, A = 0.4 and v 2 = 1.6. 

Figure 1 shows five snapshots of the wavepacket (left column) and the corresponding 
regions of the q\,p\ plane where trajectories contribute significantly to the propaga- 
tion. In these figures we have fixed the constant c = 1.0 (see equation (1401) ). except 
for figure 1(a), where c = 2.5. The width a of the smoothing Gaussian was adjusted 
to get the best results for each propagation time, starting at a = 1.5 for T = 1.0 and 
decreasing to a = 0.4 for T = 8.5 (see caption for all values). The integration over qi 
and pi was performed using a regular grid with 30 points in q 1 , varying from —3 to 
3, and 40 points in pi varying from —4 to 4. The computational time for the present 
calculation is as fast as the split-operator method, well known for being efficient and 
accurate for one dimensional problems. For T = 8.5 the calculations take about 3 
seconds in a Core 2 Quad PC with 2.4GHz. 

The wavefunctions in figure 1 were calculated using the simple discretization 



where n and m represent the grid in phase-space centered on the origin. We used a 
total of 40 and 60 points in the q and p directions respectively, with — 4 < q n < +4 
and — 6 < p m < +6. 

In spite of the accuracy of our results, specially as compared to previous calcula- 
tions using root search procedures |14j . several details remain to be understood and 
improved. The main problem is the sensitivity of the method to the choice of the 
width a and the lack of a theory on how to choose it properly and automatically. A 
possible way out of this difficulty might to be the procedure devised in [37], where 
the width is chosen to minimize the oscillations of the integrand. Another problem 
is that the propagated wavepackets turn out not to be properly normalized, and the 




(43) 



4>(B,T) = ^2( X \ Z nm)K(z* nm , Z ,T) 



AqAp 



(44) 



n,m 



amount by which normalization is lost also depends on the width a. In figure 1 the 
wave-functions have been re-normalized by hand after the propagation. 

Despite these problems the method improves the results obtained by direct compu- 
tation of the contributing trajectories and is much faster and simple to program. The 
next step is an application of the method to multidimensional systems, where the 
integrations over the initial conditions may be performed by Monte Carlo techniques. 
The difficulties just mentioned are currently under investigation. 



A Tangent matrices 



In this appendix we use the scaled units where H = b 
the tangent matrix is defined by 



1. In the u and v variables 



5u(T) 



5v{T) 



( 



\ 



M vu M vv 



( \ 

5u(0) 

v Mo) y 



(A.l) 



where 6u(0) and Sv(0) are small displacements at the initial point of the trajectory 
and 5u(T) and 5v(T) are the corresponding final deviations. The action S(v",u',T) 
for the trajectory with w(0) = u' and v (T) = v" satisfies j3] 



u(T) =u" = i 



dS 



»(0) S «- = i^. 



(A.2) 



From the differentiation of (1A.2|) keeping the variable T constant, we can obtain the 
connection between initial and final displacements. In matrix form it is 



/ \ 

5u(T) 



Sv(0) ; 



( 



\ 



q q 



( \ 

6u(0) 



(A.3) 



where S uv = d 2 S/du'dv", etc. Comparing with eq. tlA.ip we find 



.Mu 
M„ 



(A.4) 
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Using the definition of u and v in terms of q and p (notice that all these variables are 
complex) it is easy to show that [3j 

M uu = \{m qq + m pp + im pq - im qp ) 

■M uv n ijUqq Wlpp ^T^lpq ^f' '<7p) 
M vu ^(j^qq ™pp i^pq ^^qp) 
■M vv -^ijUqq TOjpp i'THpq W^qp), 

where m is the tangent matrix in the q, p system. Finally, using the definition of the 
real variables Qx, Q 2 , Pi, P 2 and defining its corresponding 4x4 tangent matrix n 
we can show that 

mq q = n n - in u 
m qp = n 13 - m i2 

(A.6) 

m pq = + in 2 i 
mpp = n 22 + in 23 . 

Therefore, by working directly with the real trajectories in the double phase space 
we can compute n and reconstruct the matrices m and M using simple linear trans- 
formations. 



B Calculation of det A 



If v(T) in equation (|29|) is an analytic function of the initial condition v%, then, by 
the Cauchy-Riemann conditions we have 

d qi (T) d Pl (T) d gi (T) d Pl (T) 

dqx dpi ' dpi dq\ 

By the definition of A, equation (|27|) . 

IF) 



On the other hand we also have, 



= i t d qi (T) + ggicm + i I'SfttT) _ a^txy 

2 y (9gi 9pi / 2 \ dpi dqi 

d gi (T) ,%(T) 
+ z- 



9gi 



9pi 



(B-3) 



and, therefore, 



det A 



dv(T) 



(B.4) 



Finally, using the second of equations 0A.20 with v' 
respect to v(T), 

d 2 S 



dvi 

dv(T) = l du'dv(T) 



V\, and differentiating with 



(B.5) 



which implies 



det A 



d 2 S 



du'dviT) 



|M„ 



(B.6) 



by equation ()A.4j) . 
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Figure Caption: 



Figure 1. (color online) The right panels show the exact and semiclassical wavefunc- 
tions for several values of T. The thin continuous line (red) displays the exact result 
obtained via split operator method; the thick solid line is the CIVR approximation 
and the dashed line (green) shows the result obtained in ref . [H] by direct computation 
of contributing trajectories. For T = 2.5 we also show the potential (V(x)/10) and 
the energy E = 2.0 of the central trajectory (shown as E/10). The left panel shows 
the contributing and non-contributing initial trajectories as white and dark areas 
respectively. The star indicates the positions q and p of the initial wavepacket. 
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